% This code can be executed to generate Figure 5A. This will plot all of
% the peak loads (dat) as a function of their junction angle. The data 
% points are colored by region and buckled conformation. In addition, 
% the average and range will be plotted in grey for all data points (datloc),
% green for nestled data points (datloc_n), and red for locked data points
% (datloc_l). The purple curve is theoretical peeling maximum (PeelF) is 
% expressed as a function of the junction angle. The yellow curve is the 
% elastic force (intersectUFelas) where the peeling and elastic energies 
% are equal (intersectU). The orange curve is point where the elastic model
% becomes linear at x=alpha (linearpts). The tensile strength of the tape 
% (TS_tape) is the black dashed line. The lap shear strength 
% (lapshearstrength) is plotted as a horizontal solid black line. The open  
% purple dots are the adjusted local peeling maxima (peelmaxdat).
% Additional formatting is done externally to match the published Firgure 5A.
% The folder titled, "Raw Data Test Files," must be open to run this code!
clear;
%Tape Constants
t=0.058; %mm
w=12.7; %mm
E=467.9745; % Pa
gamma=10.3014/1000;  %N/mm
beta=mean([189.393 190.444 184.519 182.807 186.996]);
TS_tape=23.2416;
slope=(E*w)./1000; %N/mm

%Determine Energy Equivalence and Transition to Linearity (alpha)
%initialize starting variables
angles=linspace(1.25,76,10000);
intersectUFelas=zeros(1,length(angles));
xmaxs=zeros(1,length(angles));
TS_tapeFelasx=zeros(1,length(angles));
linearpts=zeros(1,length(angles));
linearptsFelas=zeros(1,length(angles));

for i=1:length(angles)
    phi=angles(i);
    xmax=2.*w.*tand(phi./2);
    xmaxs(i)=xmax;
    %Peeling Energy
    Upeelinc=@(x) 2.*gamma.*((cscd(phi./2)).^2).*(1./2).*(x.^2).*(cotd(phi./2));
    Upeeldec=@(x) 2.*gamma.*((cscd(phi./2)).^2).*((1)./2).*((-1.*(xmax-x).^2)+2.*(xmax./2).^2).*(cotd(phi./2));
    %determine peeling energy for increaing and decreasing sections of
    %force curve
    xx=linspace(0,20,10000);
    for j=1:length(xx)
        pos=xx(j);
        if pos>=(xmax/2)
            xxinc=xx(1:j);
            Upeelincvec=Upeelinc(xxinc);
            half=j;
            break
        end
    end
    for h=half:length(xx)
        pos=xx(h);
        if pos>=xmax
            xxdec=xx(half:h);
            Upeeldecvec=Upeeldec(xxdec);
            break
        end
    end
    %add tail of constant value to make sure vectors are same length
    tailU=ones(1,(length(xx)-h)-1).*max(Upeeldecvec);
    Upeel=[Upeelincvec Upeeldecvec tailU];
    
    %Elastic Energy
    Felas=@(x) -(1/2).*E.*t.*x.*tand(phi).*log(abs((2.*x.*cotd(phi./2)-beta.*tand(phi))./(-beta.*tand(phi))));
    Felas_val=Felas(xx);
    for n=1:(length(Felas_val)-1)
        m=(Felas_val(n+1)-Felas_val(n))/(xx(n+1)-xx(n));
        b=Felas_val(n)-m*xx(n);
        %linearize elastic force
        if m>slope
            xremain=linspace(xx(n),max(xx),(length(xx)-n+1));
            Felas_val(n:length(xx))=xremain.*slope+b;
            break
        end
        linearpts(i)=xx(n);
    end
    %integrate elastic force to get elastic energy
    Uelas=cumtrapz(xx,Felas_val);
    
    %Difference between elastic and peeling energies to determine where
    %they become equal
    subU=Uelas-Upeel;
    for k=1:length(subU)-1
        test=subU(k);
        if test>0
            intersectUFelas(i)=Felas_val(k);
            break
        end
    end

    %Determine Alpha
    linearptsFelas(i)=Felas_val(n);
end
anglestheo=angles;

%Interperet Raw Data: different sets of data require different treatment
%due to the organization in the .txt files and category of the sample that
%the file cooresponds to.

%Tensile Tester Trials
names=['2019-02-08_001_70.0.txt',...
    '2019-02-27_002_75.0.txt',...
    '2019-02-27_003_66.0.txt',...
    '2019-02-27_004_65.0.txt',...
    '2019-02-27_005_52.6.txt',...
    '2019-02-27_006_40.9.txt',...
    '2019-02-27_007_28.0.txt',...
    '2019-03-01_008_85.1.txt',...
    '2019-03-01_009_52.0.txt',...
    '2019-03-01_010_49.5.txt',...
    '2019-03-01_011_30.6.txt',...
    '2020-04-15_001_10.0.txt',...
    '2020-04-15_002_12.5.txt',...
    '2020-04-15_003_15.0.txt',...
    '2020-04-15_004_17.5.txt',...
    '2020-04-15_005_20.0.txt',...
    '2020-04-15_006_22.5.txt',...
    '2020-04-15_007_25.0.txt',...
    '2020-04-15_008_27.5.txt',...
    '2020-04-15_009_30.0.txt',...
    '2020-04-15_010_32.5.txt',...
    '2020-04-15_011_35.0.txt',...
    '2020-04-15_012_37.5.txt',...
    '2020-04-15_013_40.0.txt',...
    '2020-04-15_014_42.5.txt',...
    '2020-04-15_015_45.0.txt',...
    '2020-04-15_016_47.5.txt',...
    '2020-04-15_017_50.0.txt',...
    '2020-04-15_018_52.5.txt',...
    '2020-04-15_019_55.0.txt',...
    '2020-04-15_020_57.5.txt',...
    '2020-04-15_021_60.0.txt',...
    '2020-04-15_022_65.0.txt',...
    '2020-04-15_023_70.0.txt',...
    '2020-04-15_024_75.0.txt',...
    '2020-05-04_025_10.0.txt',...
    '2020-05-04_026_12.5.txt',...
    '2020-05-04_027_15.0.txt',...
    '2020-05-04_028_17.5.txt',...
    '2020-05-04_029_20.0.txt',...
    '2020-05-04_030_22.5.txt',...
    '2020-05-04_031_25.0.txt',...
    '2020-05-04_032_27.5.txt',...
    '2020-05-04_033_30.0.txt',...
    '2020-05-04_034_17.5.txt',...
    '2020-05-04_035_20.0.txt',...
    '2020-05-04_036_22.5.txt',...
    '2020-05-04_037_25.0.txt',...
    '2020-05-04_038_27.5.txt',...
    '2020-05-04_039_30.0.txt',...
    '2020-05-04_040_20.0.txt',...
    '2020-05-04_041_22.5.txt',...
    '2020-05-04_042_25.0.txt',...
    '2020-05-04_043_27.5.txt',... 
    '2020-05-12_046_17.5.txt',...
    '2020-05-12_047_20.0.txt',...
    '2020-05-12_048_22.5.txt',...
    '2020-05-12_049_25.0.txt',...
    '2020-05-12_050_27.5.txt',...
    '2020-05-12_051_30.0.txt',...
    '2020-05-12_052_17.5.txt',...
    '2020-05-12_053_20.0.txt',...
    '2020-05-12_054_22.5.txt',...
    '2020-05-12_055_25.0.txt',...
    '2020-05-12_056_27.5.txt',...
    '2020-05-12_057_30.0.txt',...
    '2020-05-12_058_17.5.txt',...
    '2020-05-12_059_20.0.txt',...
    '2020-05-12_060_22.5.txt',...
    '2020-05-12_061_25.0.txt',...
    '2020-05-12_062_27.5.txt',...
    '2020-05-12_063_30.0.txt',...
    '2020-05-12_064_17.5.txt',...
    '2020-05-12_065_20.0.txt',...
    '2020-05-12_066_22.5.txt',...
    '2020-05-12_067_25.0.txt',...
    '2020-05-12_068_27.5.txt',...
    '2020-05-12_069_30.0.txt',...
    '2020-05-14_070_10.0.txt',...
    '2020-05-14_071_12.5.txt',...
    '2020-05-14_072_15.0.txt',...
    '2020-05-14_073_10.0.txt',...
    '2020-05-14_074_12.5.txt',...
    '2020-05-14_075_15.0.txt',...
    '2020-05-14_076_17.5.txt',...
    '2020-05-21_077_15.0.txt',...
    '2020-05-21_078_42.5.txt',...
    '2020-05-21_079_35.0.txt',...
    '2020-05-21_080_50.0.txt',...
    '2020-05-21_081_37.5.txt',...
    '2020-05-21_082_45.0.txt',...
    '2020-05-21_083_47.5.txt',...
    '2020-05-21_084_32.5.txt',...
    '2020-05-21_085_40.0.txt',...
    '2020-05-21_086_22.5.txt',...
    '2021-03-16_100_50.0.txt',...
    '2021-03-16_101_15.0.txt',...
    '2021-03-16_102_15.0.txt',...
    '2021-03-16_103_15.0.txt',...
    '2021-03-16_104_17.5.txt',...
    '2021-03-16_105_10.0.txt'];
%short algorithm to determine where each file name starts and ends to
%determine the junction angle information which is stored in the filenames
for i=1:length(names)
    chunk=i+3;
    test=names(i:chunk) == '.txt';
    if all(test)
        numchar=chunk;  % number of characters in each file name
        break
    end
end
%initialize variables for data analysis
ntrials=length(names)/numchar;
jj=zeros(1,ntrials);
jj(1)=1;
kk=zeros(1,ntrials);
kk(1)=numchar;
elasticmax=[];
elasticmaxangles=[];
peelmax=[];
peelmaxangles=[];
peakload=[];
angles=[];
for i=1:ntrials
    clear Data displacement load time
    filename=names(jj(i):kk(i));
    fileID=importdata(filename);
    j=jj(i)+numchar;
    k=kk(i)+numchar;
    jj(i+1)=j;
    kk(i+1)=k;
    Data=fileID.data;
    if str2num(filename(4))==9  %older data sets require different treatment
        displacement=Data(:,1).*(50E-3);  %mm
        load=abs(Data(:,2)-abs(max(Data(:,2))));  %N
        adjusttrials=[2 3 5 6 10];       %we need to pick out the elastic and peeling peaks.
        adjusttripletroubletrials=[10];  %theoretically the highest peak is the elastic peak and the second is the peeling
        adjustquadtroubletrials=[];      %however, some trials have multiple peaks in the data so we basically manually pickout
        adjustbigtroubletrials=[2 3];    %the proper peaks. Only 5 samples have this issue and are labeled by the 
    else                                 %"adjust#troubletrials" varialbes
        displacement=Data(:,1);  %mm
        load=Data(:,2);  %N
        adjusttrials=[13 15 17 18 19 20 22 24 81 82 85];
        adjusttripletroubletrials=[81];
        adjustquadtroubletrials=[19];
        adjustbigtroubletrials=[];
    end
    phi=str2num(filename(16:19));
    %begin to pick out peaks
    if any(ismember(adjusttrials,str2num(filename(12:14)))) %these trials have 2 peaks, global: elastic, local: peeling
        [pks, pkloc]=findpeaks(load,displacement);
        [pks2, pkloc2]=findpeaks(pks,pkloc);
        allpks=[pks2 pkloc2];
        sortpks=sortrows(allpks,'descend');
        peelmax=[peelmax sortpks(2,1)];
        peelmaxangles=[peelmaxangles phi];
        elasticmax=[elasticmax sortpks(1,1)];
        elasticmaxangles=[elasticmaxangles phi];
        if any(ismember(adjusttripletroubletrials,str2num(filename(12:14)))) 
            peelmax(length(peelmax))=sortpks(3,1);  %these trials have the issues discussed above
        end
        if any(ismember(adjustquadtroubletrials,str2num(filename(12:14))))
            peelmax(length(peelmax))=sortpks(4,1);
        end
        if any(ismember(adjustbigtroubletrials,str2num(filename(12:14))))
            [pks3, pkloc3]=findpeaks(pks2,pkloc2);
            allpks=[pks3 pkloc3];
            sortpks=sortrows(allpks,'descend');
            if str2num(filename(12:14))==2
                peelmax(length(peelmax))=sortpks(2,1);
            end
            if str2num(filename(12:14))==3
                peelmax(length(peelmax))=sortpks(3,1);
            end
        end
    else
        [pks, pkloc]=findpeaks(load,displacement);  %these trials are peeling only
        allpks=[pks pkloc];
        sortpks=sortrows(allpks,'descend');
        peakload=[peakload sortpks(1,1)];
        angles=[angles phi];
    end
end

%Peak Load Only Trials
peakloaddata=importdata('Peak Load Only Raw Data.txt');
peakloaddata2=peakloaddata.data;
anglesexpinv=peakloaddata2(:,1);
maxloadexp=flip(peakloaddata2(:,2))';
anglesexp=flip(anglesexpinv)';

% Analysis

%organize normalized data into single matrix
dat=[anglesexp angles; maxloadexp./areacalculation(w,anglesexp) peakload./areacalculation(w,angles)];
dat=sortrows(dat',1);
elasticmaxdat=[elasticmaxangles; elasticmax./areacalculation(w,elasticmaxangles)];
peelmaxdat=[peelmaxangles; peelmax./areacalculation(w,peelmaxangles)];

%Peeling Model
PeelF=@(phi) 2.*gamma.*w.*((cscd(phi./2)).^2)+0.4536;
phiphi=linspace(min(dat(:,1)),90,10000);
peel_f=PeelF(anglestheo)./areacalculation(w,anglestheo);

%Determine Boundary Between Regions II and III and Regions I and II
linearptsFelas=linearptsFelas./areacalculation(w,anglestheo);
intersectUFelas=intersectUFelas./areacalculation(w,anglestheo);
for i=1:length(anglestheo)
    alphaloc=linearptsFelas(i);
    peelloc=peel_f(i);
    if peelloc<alphaloc
        R_I_II=(anglestheo(i)+anglestheo(i-1))/2;
        break
    end
end
for i=1:length(anglestheo)
    peelloc=peel_f(i);
    eqUloc=intersectUFelas(i);
    if anglestheo(i)>R_I_II && peelloc>eqUloc
        R_II_III=(anglestheo(i)+anglestheo(i-1))/2;
        break
    end
end

%Plotting
figure;
hold on

%Determine Regions & Plot Data Points
R1loads=[];
R2loads_l=[];
R2loads_n=[];
R3loads_l=[];
R3loads_n=[];
R1angles=[];
R2angles_l=[];
R2angles_n=[];
R3angles_l=[];
R3angles_n=[];

for c=1:length(dat(:,1))
    phi=dat(c,1);
    if phi<21.9
        scatter(dat(c,1),dat(c,2),30,'filled','MarkerFaceAlpha',1/3,'MarkerFaceColor',[0 0 0]...
            ,'MarkerEdgeColor',[0 0 0],'MarkerEdgeAlpha',2/3);
        R1loads=[R1loads dat(c,2)];
        R1angles=[R1angles dat(c,1)];
    elseif phi<=31 && dat(c,2)>=0.0256
        scatter(dat(c,1),dat(c,2),30,'filled','MarkerFaceAlpha',1/3,'MarkerFaceColor',[255/255 112/255 101/255]...
            ,'MarkerEdgeColor',[0 0 0],'MarkerEdgeAlpha',2/3);
        R2loads_l=[R2loads_l dat(c,2)];
        R2angles_l=[R2angles_l dat(c,1)];
    elseif phi<=31 && dat(c,2)<0.0256
        scatter(dat(c,1),dat(c,2),30,'filled','MarkerFaceAlpha',1/3,'MarkerFaceColor',[0/255 190/255 30/255]...
            ,'MarkerEdgeColor',[0 0 0],'MarkerEdgeAlpha',2/3);
        R2loads_n=[R2loads_n dat(c,2)];
        R2angles_n=[R2angles_n dat(c,1)];
    else
        scatter(dat(c,1),dat(c,2),30,'filled','MarkerFaceAlpha',1/3,'MarkerFaceColor',[0/255 190/255 30/255]...
            ,'MarkerEdgeColor',[0/255 100/255 0/255],'MarkerEdgeAlpha',2/3);
        R3loads_n=[R3loads_n dat(c,2)];
        R3angles_n=[R3angles_n dat(c,1)];
    end
end
scatter(elasticmaxdat(1,:)',elasticmaxdat(2,:)',30,'filled','MarkerFaceAlpha',128/255,'MarkerFaceColor',[255/255 112/255 101/255]...
    ,'MarkerEdgeColor',[180/255 0/255 0/255],'MarkerEdgeAlpha',128/255);
R3loads_l=[R3loads_l elasticmaxdat(2,:)];
R3angles_l=[R3angles_l elasticmaxdat(1,:)];

scatter(peelmaxdat(1,:)',peelmaxdat(2,:)',30,'filled','MarkerFaceAlpha',0/3,'MarkerFaceColor',[180/255 0/255 0/255]...
    ,'MarkerEdgeColor',[127/255 0/255 255/255],'MarkerEdgeAlpha',3/3);
R3loads_n=[R3loads_n peelmaxdat(2,:)];
R3angles_n=[R3angles_n peelmaxdat(1,:)];

nlocked=length(R3loads_l)+length(R2loads_l);
nnestled=length(R3loads_n)+length(R2loads_n);
percentlocked=(nlocked/nnestled);

%Determine Mean and Range for given sections of the peeled data to draw the 
%Mean and Range lines and regions respectively
rangedat=sortrows([dat(:,1)' elasticmaxdat(1,:); dat(:,2)' elasticmaxdat(2,:)]',1);

binlims=[linspace(min(rangedat(:,1)),20,5) linspace(22,30,4) linspace(33,max(rangedat(:,1)),9)];
numbin=length(binlims);
datloc=zeros(6,numbin-1);
for i=1:(numbin-1)
    minlim=binlims(i);
    maxlim=binlims(i+1);
    if i==1
        datloc(1,i)=minlim;
    elseif i==(numbin-1)
        datloc(1,i)=maxlim;
    else
        datloc(1,i)=(maxlim+minlim)/2;
    end
    bin=[];
    for j=1:length(rangedat(:,1))
        angle=rangedat(j,1);
        if angle>=minlim && angle<=maxlim
            bin=[bin rangedat(j,2)];
        end
    end
    datloc(2,i)=mean(bin);
    datloc(3,i)=std(bin);
    if length(bin)==1
        datloc(4,i)=mean(bin);
        datloc(5,i)=mean(bin);
    else
        datloc(4,i)=min(bin);
        datloc(5,i)=max(bin);
    end
    datloc(6,i)=length(bin);
end

%Plot Range and Average for All Data Points
jbfill(datloc(1,:),datloc(4,:),datloc(5,:),[200/255 200/255 200/255],'none',[0 0 0],119/255);
plot(datloc(1,:),datloc(2,:),'Marker','none','Color',[130/255 130/255 130/255],'LineStyle','-','LineWidth',2);

%optional: plot lines to mark the range
%plot(datloc(1,:),datloc(4,:),'Marker','none','Color',[130/255 130/255 130/255],'LineStyle','-');
%plot(datloc(1,:),datloc(5,:),'Marker','none','Color',[130/255 130/255 130/255],'LineStyle','-');


%Determine Range and Average for Nestled Data Points
rangedat_n=sortrows([R2angles_n R3angles_n; R2loads_n R3loads_n]',1);

binlims_n=[linspace(20,30,3) linspace(33,75,6)];
numbin_n=length(binlims_n);
datloc_n=zeros(6,numbin_n-1);
for i=1:(numbin_n-1)
    minlim=binlims_n(i);
    maxlim=binlims_n(i+1);
    if i==1
        datloc_n(1,i)=minlim;
    elseif i==(numbin_n-1)
        datloc_n(1,i)=maxlim;
    else
        datloc_n(1,i)=(maxlim+minlim)/2;
    end
    bin=[];
    for j=1:length(rangedat_n(:,1))
        angle=rangedat_n(j,1);
        if angle>=minlim && angle<=maxlim
            bin=[bin rangedat_n(j,2)];
        end
    end
    datloc_n(2,i)=mean(bin);
    datloc_n(3,i)=std(bin);
    if length(bin)==1
        datloc_n(4,i)=mean(bin);
        datloc_n(5,i)=mean(bin);
    else
        datloc_n(4,i)=min(bin);
        datloc_n(5,i)=max(bin);
    end
    datloc_n(6,i)=length(bin);
end

%Plot Range and Average for Nestled Data Points
jbfill(datloc_n(1,:),datloc_n(4,:),datloc_n(5,:),[0/255 190/255 30/255],[1 1 1],[0 0 0],50/255);
plot(datloc_n(1,:),datloc_n(2,:),'Marker','none','Color',[0/255 100/255 0/255],'LineStyle','-','LineWidth',2);

%optional: plot lines to mark the range
%plot(datloc_n(1,:),datloc_n(4,:),'Marker','none','Color',[0/255 100/255 0/255],'LineStyle','-');
%plot(datloc_n(1,:),datloc_n(5,:),'Marker','none','Color',[0/255 100/255 0/255],'LineStyle','-');


%Determine Range and Average for Locked Data Points
rangedat_l=sortrows([R2angles_l R3angles_l; R2loads_l R3loads_l]',1);

binlims_l=[linspace(20,30,3) linspace(33,75,6)];
numbin_l=length(binlims_l);
datloc_l=zeros(6,numbin_l-1);
for i=1:(numbin_l-1)
    minlim=binlims_l(i);
    maxlim=binlims_l(i+1);
    if i==1
        datloc_l(1,i)=minlim;
    elseif i==(numbin_l-1)
        datloc_l(1,i)=maxlim;
    else
        datloc_l(1,i)=(maxlim+minlim)/2;
    end
    bin=[];
    for j=1:length(rangedat_l(:,1))
        angle=rangedat_l(j,1);
        if angle>=minlim && angle<=maxlim
            bin=[bin rangedat_l(j,2)];
        end
    end
    datloc_l(2,i)=mean(bin);
    datloc_l(3,i)=std(bin);
    if length(bin)==1
        datloc_l(4,i)=mean(bin);
        datloc_l(5,i)=mean(bin);
    else
        datloc_l(4,i)=min(bin);%mean(bin)-std(bin);% can also plot std instead of range
        datloc_l(5,i)=max(bin);%mean(bin)+std(bin);%
    end
    datloc_l(6,i)=length(bin);
end

%Plot Range and Average for Locked Data Points
jbfill(datloc_l(1,:),datloc_l(4,:),datloc_l(5,:),[255/255 112/255 101/255],[1 1 1],[0 0 0],50/255);
plot(datloc_l(1,:),datloc_l(2,:),'Marker','none','Color',[180/255 0/255 0/255],'LineStyle','-','LineWidth',2);

%optional: plot lines to mark the rang
%plot(datloc_l(1,:),datloc_l(4,:),'Marker','none','Color',[180/255 0/255 0/255],'LineStyle','-');
%plot(datloc_l(1,:),datloc_l(5,:),'Marker','none','Color',[180/255 0/255 0/255],'LineStyle','-');


%Plot Peeling Model
plot(anglestheo,peel_f,'Marker','none','Color',[127/255 0 255/255],'LineStyle','-','LineWidth',2);

%Plot Tensile Strength of Tape
hline=ones(1,length(phiphi)).*TS_tape;
plot(phiphi,hline./areacalculation(w,phiphi),'Color','k','LineStyle','--');


%Determine & Plot Lap Shear Strength
LSSdata=importdata('LSS Raw Data.txt');
Data=LSSdata.data;
JL=Data(:,1)';
PeakForce=Data(:,2)';
JunctionArea=(JL).*w;
ff=@(m,x) m.*x;
[linreg]=fit(JunctionArea',PeakForce','poly1');
lapshearstrength=linreg.p1;
LSF=ones(1,length(phiphi)).*lapshearstrength;
plot(phiphi,LSF,'k');

%Plot where Elastic Energy=Peeling Energy
plot(anglestheo,intersectUFelas,'Color',[255/255 193/255 37/255],'LineWidth',2);

%Plot where slope of elastic force=modulus of single tape (alpha)
plot(anglestheo,linearptsFelas,'Color',[255/255 126/255 24/255],'LineWidth',2);

%Plot Region Boundaries
line1_2=ones(1,1000).*R_I_II;
line2_3=ones(1,1000).*R_II_III;
yy=linspace(0,0.1,1000);
plot(line1_2,yy,'b','LineStyle','--');
plot(line2_3,yy,'b','LineStyle','--');

%Formatting the plot
xlabel('Junction Angle (\circ)','FontSize',14); 
xlim([10 90]);
ylabel('Max Load/Junction Area (MPa)','FontSize',14);
ylim([0 0.1]);
ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XAxis.Limits(1):2.5:ax.XAxis.Limits(2);
ax.FontSize=14;